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NUMERICAL  SOLUTION  OF  A  SYSTEM  OF  LANCH ESTER 
DIFFERENTIAL  EQUATIONS 

by 

Svein  B.  Haugberg 


ABSTRACT 


Using  the  Runge-Kutta-Merson  method  a  computer  program  was 
developed  for  the  numerical  solution  of  a  system  of  generalized 
Lanchester  equations.  A  detailed  description  is  given  of  the 
input-output  features  of  the  computer  program  and  of  the  algorithm 
used  to  solve  the  equations.  Two  sample  problems  illustrate  the 
use  of  the  program  and  its  flexibility. 


1 


INTRODUCTION 


The  purpose  of  this  memorandum  is  to  describe  the  development  of  a 
computer  program  for  numerically  solving  a  system  of  generalized 
Lanchester  equations. 

The  system  of  generalized  Lanchester  differential  equations  [Ref.  1] 
may  be  written  as  follows; 

X  »  -f (X)  +  r  ,  [Eq.  1] 

where  f_(X)  ^0  for  X  so,  and  with  known  initial  values  X(0). 

The  vector  X  ”  fXj,  X2  ,  .  Xn) ,  X^  -  X^(t),  describes  the  number 
of  units  engaged  in  combat  in  each  of  n  different  classes  as  a 
function  of  time.  The  components  tVCX)  of  £(X)  are  the 
expected  numbe*-  of  units  of  class  i  lost  per  unit  time.  The 
vector  r  -  ( r;  ,  rS)...,  r^) ,  r,  =r.  (t;),  is  the  rate  of 
introduction  of  new  units  in  each  of  the  n  classes.  Each 
f|QC)  is  assumed  to  be  the  sum  of  nu  contributions; 

mi 

MX)  *  T.  r.  .  (X)  ,  [Eq.  2] 

j,=l  u 

i,  L 

where  i\  .(X)  Is  the  j  '  contribution  to  £\(X).  These 
contributions  are  generally  associated  with  different  factors 
causing  attrition  of  type  i  units. 

The  remainder  of  this  memorandum  proceeds  along  the  following 
lines.  First,  the  method  of  numerically  solving  the  system  of 
Lanchester  equations,  called  the  Runge-Kutta-Merson  Method,  is 
discussed.  Second,  a  computer  program  that  performs  the  numerical 
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solution  is  described.  Finally,  two  sample  problems  and  their 
computer  solutions  are  presented.  A  complete  listing  of  the 
computer  program  appears  in  Appendix  A. 


1. 


THE  RUNGE-Kl' TTA-MERSON  METHOD  CF  SOLUTION 


The  Runge-Kutta  method  of  numerical  integration  is  nne  of  several 
fourth-order  processes  and  is  applicable  to  all  systems  of  1st 
order  differential  equations  of  the  form  X(t)  -  g^[X(t),  t],  and 
therefore  also  to  all  systems  of  higher  order  that  can  be  reduced 
to  the  same  form.  To  start  the  integration,  initial  values  cf 
t  o )  X(t0)  are  required.  The  method  employs  a  constant  step  length 
in  the  integration.  A  difficulty  with  this  method  is  that  it  is 
often  not  clear  how  to  determine  a  proper  step  size  to  achieve  a 
known  local  truncation  error  at  each  step.  The  method  of  Mer-son 
overcomes  this  difficulty  by  automatically  determining  the  step 
length  required  to  obtain  a  predetermined  accuracy, 

If  we  define  £(X,  t)  to  be 

£(X,  t)  ---  -fI  X)  *•  r  ,  [Eq.  1] 

then  the  system  of  I.anehester  equations  described  above  becomes 


X  -■*  £(X,  t)  .  [Eq.  4  ] 

This  more  compact  form  will  be  used  for  the  following  discussion 
of  the  Runge-Kutta-Merson  method . 

Given  t0,  an  initial  time,  >  the  corresponding  initial 

value  of  X,  and  h,  the  step  length,  the  Runge-Kutta-Merson  method 


[Ref ,  2]  allows  one  to 

i ntegratc 

g(X,  t)  to 

a  Later 

t  i  me 

tj  =  t0  +  h  and  obtain 

X(  t,  ) 

hy 

the  following 

ptoredu 

re  . 

First,  let  t.|t  t  ^ ,  t( 

and 

Si 

bo  four  inter 

■meri  i  at  e 

points  in 

one  integration  step  h  from  t0  to  tj  selected  as  follows: 
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Th<*  inothod  of  Mrison  lot  .tut.om.it  i  <■;»  I  I  y  adjusting  tin*  step  I onyt h  h 
is  as  follows. 


Given  a  predetermined  relative  accuracy  e,  if  q^  >  e  for  at 
least  one  i,  the  step  length  will  be  halved  and  the  integration 
repeated  with  this  new  h.  This  will  go  on  until  <  e  for  all  i. 
t f*  ^  on  th?  of Ku n  h 2nd f  c 1  _  «*"  «  /  o  o  fcr*  2 1 1  i.  2 vC p  1  ci* 'Tth  v» i  1 1 

be  doubled,  and  this  new  h  will  be  used  in  the  next  step. 

To  continue  the  integration,  let  tj  and  X^l  )  be  the  new 
initial  values  t0  and  X(t0)  respectively,  and  repeat  the 
procedure  described  above. 

The  effect  of  non-linearities  in  g  may  overestimate  the  smallness 
of  the  step  length  necessary  to  achieve  the  required  accuracy. 

It  is  also  possible  that  the  Merson  process  might  underestimate 
the  error,  but  according  to  Ref.  2  no  example  of  this  type  of 
behaviour  has  been  encountered. 

Applied  to  a  computer  program,  the  advantage  of  the  Runge-Kutta-Merson 
method  consequently  is  that  the  integration  always  will  be  performed 
with  the  predetermined  accuracy  and  that  the  step  length  is  never 
much  smaller  than  necessary,  which  saves  computer  time. 
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2. 


COMPUTER  PROGRAM 


2. 1  Introduction 


A  computer  program  that  performs  numerical  solutions  of  the  system 
of  generalized  Lanchester  equations  described  in  the  previous 
section  has  been  written  in  ALGOL  for  the  Elliott  503  computer. 

The  results  produced  by  the  program  are  the  values  of  X  and 
the  other  time-dependent  variables  in  Eqs.  1  and  2, 

for  certain  discrete  values  of  time  in  a  specified  region  [t^n>  ^m]j 
where  t^n  is  the  initial  time ,  X(t^n)  is  assumed  known, 

The  program  is  built  around  a  procedure  MERSON,  based  on  the  Runge- 
Kutta-Merson  method.  For  given  values  of  X  at  a  time  t,  the 
procedure  MERSON  computes  £  at  a  later  time  t  +  Dt  by 
successively  increasing  time  with  variable  step  length  h  until 
t  +  Dt  is  reached. 

Another  procedure,  FUNC,  contains  all  the  functions  fjj(X)  and 
r^(t),  and  this  procedure  computes  the  values  of  the  derivatives  X 
when  called  by  MERSON. 

2. 2  The  Main  Program 

The  central  part  of  the  main  program  consists  of  a  loop  that 

starts  with  the  initial  values  t.  and  X(t.  )  and  computes  X 

in  —  an  — 

at  a  later  time  t  = t^n  +  Dt  by  calling  MERSON.  The  values  of  X 

at  this  new  time  arc  used  to  compute  the  other  time-dependent 

variables  for  this  t.  These  are  obtained  by  calling  FUNC.  The 

values  of  all  the  time-dependent  variables  are  stored  in  arrays 

for  later  output.  In  the  next  run  through  the  loop  the  last  values 

of  t  and  X  are  new  initial  values,  and  MERSON  computes  X  at 

t  +  Dt.  The  loop  terminates  at  t~t  (t  -  t.  is  an  integer 

m  m  in  b 

mu  It ipl e  of  Dt ) . 


The  rest  of  the  main  program  provides  for  output.  Some  flexibility 
exists  in  the  choice  of  variables  output;  the  different 
possibilities  will  be  given  in  the  description  of  input  and 
output . 

2 . 3  The  Procedures  MERSON  and  FUNC 


The  Runge-Kutta-Merson  method  is  described  in  Chapter  1.  This 
method  is  well  suited  for  implementation  on  an  electronic  computer 
[Ref.  3].  Equations  5,  6  and  7,  together  with  the  accuracy 

test, const itute  the  main  part  of  the  procedure  MERSON.  The  following 
quantities  are  required  by  the  procedure  MERSON  when  it  is  called 
in  the  main  program  to  compute  X(t;  +  Dt)  for  known  X(t)! 

n,  number  of  equations. 

t,  start  point  of  integration. 

-  t  +  Dt,  end  point  of  integration. 

x.(t),  (  i  =  1 ,  n )  ,  the  values  of  X  in  the  start  point 

1  of  integration. 

h,  a  starting  value  for  the  step  length. 

hmin,  lower  bound  for  the  step  length.  (This  prevents 

indefinite  cycling,  otherwise  possible  in  some  cases.) 

e,  desired  bound  for  relative  accuracy. 

With  these  quantities  given, the  procedure  MERSON  performs  one  step  h 
according  to  Eqs.  b  and  7  and  the  accuracy  test  with 
t0  ‘  t,  and  calculates  X(t}  )  .  The  derivatives  X.  ~  g.  (i  1,  n)  , 
needed  in' every  substep,  are  calculated  by  calling  FliNC  before  each’ 
of  them.  The  computations  are  continued  by  letting  tj  and  )£(t.j  ) 
be  tin:  new  t0  and  X(t0)  respectively.  When  the  integration  has 
reached  a  stage  where  t0+  h  exceeds  the  end  point  of  integration, 

/,  a  final  step  length  /-  t0  will  be  used  to  reach  /  exactly. 

For  this  last,  special  step  a  better  accuracy  than  the  one  needed 
will  result  because  s  -  tQ<h. 


Compute  for  the  5th  substep,  and 
the  corresponding  truncation  error 


for  all  i 
from  1  to  n 


Is  the 

"truncation  error  greatei 
than  the  desired 
'"'~"~-^ac  c  u  r  a  c  y 


If  h  is  less  than  or  equal  to  h  min, 
print  the  relative  truncation  error 


t  is  now  equal  to  the  time  at  the 
end  point  of  the  step  h 


Halve  the  step  length 


"Is  h  less  ^ 
than  h  min?. 


the  truncations^ 
error  less  than  1/32  of 
^  the  desi red  accuracy 
for  all  i  ? 


Set  h  equal 
to  h  min 


Restore  t  and  X(t)  and  go  back 
for  repeated  integration  with 
this  new  step  length 


1ft  is  not  equal  No' 
to  the  end  point  z 
double  the  step 
length 


■~^Is  t  equal  to 
^hi e  end  point  z  ? 


2—@ 


End  MERSON 


FiG  1  FLOW  Di  AGRAM  FOR  THE  PROCEDURE  MERSON 


Results  produced  by  MERSON  are: 


t,  its  value  is  now  that  ot'  /  ~  t  f  Dt . 

X.(t),  ti  =  1,  n),  the  values  of  X  at  the  end  point. 

^  « 
n  .  (  x .  t).  (i—  1.  n  I  .  the  values  of  the  H*>  e  i  v  at  i  vo  *  X  ^ 

jL  at  the  end  point, 
h,  adjusted  step  length. 

A  flow  diagram  for  the  procedure  MERSON  is  shown  in  Fig.  1. 

It  appears  from  the  flow  diagram  that  if  h  becomes  less  than 
hmin,  h  min  is  used  as  the  step  length  instead  of  h.  This  is  in 
order  to  avoid  particularly  long  run  time  on  the  computer.  In  this 
case,  however,  it  might  happen  that  the  desired  accuracy  is  not 
satisfied.  As  an  indication  for  the  user  of  the  program,  the 
character  followed  by  the  relative  truncation  error  is  printed 
in  the  table  of  X^(t)  every  time  hmin  is  used  as  the  step 
1 ength . 

1 . 4  Output  and  Input. 

The  results  of  the  computations  are  printed  out  for  certain  discrete 

values  of  t,  starting  with  tjn  and  increasing  repeatedly  with  a 

constant  time  interval  Dt  until  t  is  reached.  Since  the 

m 

integration  of  each  interval  Dt  is  an  independent  computation 
with  new  initial  conditions,  these  initial  conditions  can  be  changed 
at  the  start  point  of  any  interval.  Therefore  at  these  points  the 
analytical  form  of  the  functions  f .  .  (_X)  and  r.(t-)  ran  be  changed, 
and  instantaneous  i nt ruduct i ons  of  units  of  any  class  i  can  be  made 
Whenever1  such  changes  or  introductions  occur,  the  region  of 
integration  has  to  he  divided  into  smaller  regions  such  that  all 
change's  and  introductions  occur  at  the  start  point,  of  a  region,  and 
that,  every  region  is  an  integral  number  Dt . 


0  ut,  put 

As  mentioned  above,  the  program  allows  for  a  choice  in  t he  selection 
of  variables  and  combinations  of  them  given  as  output.  The  output 


In 


cons i st s  of : 


a.  The  number  of  units  introduced  instantaneously  and  the 
corresponding  instants  in  time. 


b. 

Xi 

VS 

t  for  i  —  1 ,  n  . 

c . 

X. 

1 

vs 

Xj  for  specified  index 

pairs  (i,  j). 

d. 

f  . 

1 

vs 

t  for  specified 

indices 

i  . 

e . 

r . 

i 

V.S 

t  for  specified 

indices 

i  . 

f . 

f.  ./f. 

(in  percent)  v-s 

t  for 

specified  indices 

and  specified  indices  j  for  each  i. 
g.  f^/fj  vs  t  for  specified  index  pairs  (i,  j). 


The  last  five  outputs  are  optional  and  may  be  skipped  by  not 
specifying  any  indices.  This  output  form  enables  the  user  of  the 
program  to  study  the  importance  of  instantaneous  introductions  and 
changes  in  attrition  functions  and  introduction  rates,  and  the 
relative  importance  of  the  different  contributions  f ^ ^ ( X)  causing 
attrition  of  type  i  units.  This  will  be  illustrated  by  the 
sample  problems  in  the  next  chapter. 

There  are  two  forms  of  output:  paper  tape  and  plotted.  The  paper 
tape  output  is  punched  on  the  Number  1  punch,  and,  of  course,  must 
be  run  through  the  Flexowriter  to  obtain  a  printer  output.  The 
plotted  output  is  generated  by  the  Caleomp  plotter. 

Input 

There  are  two  kinds  of  input  to  this  program:  input  parameters 
and  input  functions.  These  are  listed  in  Tables  J  and  2. 


1  i 


TABLE  1 

LIST  OF  INPUT  PARAMETERS 


In  the  following  list  the  letter-  in  Mrcnrhp««  following  the 
explanation  indicates  whether  the  variable  is  integer  (I)  or 
real  (R) . 


N,  number  of  differential  equations  (1) 

DT,  the  constant  time  interval  between  two  following  (R) 

output  points 

II ,  a  start  value  for  the  step  length  (R) 

II  MIN,  lower  bound  for  the  step  length  (R) 

F,  desired  bound  for  relative  accuracy  (R) 

MI,  maximum  number  of  contributions  of 

f  .  .  (X)  to  f  .  (  X)  ( i  l ,  n ) ,  [  M I  "  max  (m^  ) »  i  “  1 ,  n  ]. 

N I ,  number  of  time  regions  according  to  the  instantaneous  (I) 

introductions  and  the  change  of  functions 

NS,  number  of  pairs  ( X ^ ,  X^)  to  be  plotted  (I) 

(  X  ^  vs  x.,  t  as  a  parameter) 

Nl,  number  of  outputs  (print  and  plot)  of  f ^  Vs  t 

NR,  number  of  outputs  (plot)  of  r.  Vs  t  (1) 

N2,  number  of  outputs  (print  and  plot)  of 

f,  ,/f.  vs  t  (number  of  different  indices  i) 

'  .1  ' 

NP,  number  of  out-puts  (print  and  plot)  of  (1) 

MA[  L  (L  I,  Nil,  number  of  Dt-  in  each  of  NI  regions  (l) 

(the  lengths  of  the  regions  are  MAjl.  DT). 

I T  ,  initial  time 

DXfl.,,1]  (I  I,  N  for  each  I  1.  Nil.  i  nst  ant  aneous  IR) 

i  nt.  roduet  i  oils  til  units  in  the  start  point  of 
region  l. .  The  initial  value*,  of  Xjti  1.  nl 
are  i  nr  1  udiul  here  (1  II, 


I' ,  V  (K  -  1,  NS;,  specification  of  indices  in  NS  plots  (11 

of  X  .  vs  X  . 
i  J 

f  (1=1,  Nl),  specification  of  index  in  Ni  outputs  (I) 

Ul  I  j  V-»  1 

U,  Ml  (J=  1,  N2),  specification  of  index  i  and  number  of  (II 

different  indicox  j  for  each  i  in  outputs 

of  f .  .  /  f  .  v->  t 

i  J  i 

V  (K  =  l,  Ml  for  each  J  =  1 ,  N2),  specification  of  index  j  (1) 

in  Ml  outputs  of  f.j/f.  vs  t 

U,  V  ( J  —  J. ,  NP),  specification  of  indices  i  and  j  in  NP  (I) 
outputs  of  t\/f  .  vs  t, 

I  (II,  NR),  specification  of  index  in  NR  outputs  (I) 

o  t  r .  vs  t 

i 

The  input  data  must  be  given  numerical  values  and  punched  on  a 
data  tape  in  the  following  order: 

N,  DT,  H,  H  MIN,  E,MI,  NI,  NS,  Nl,  NR,  N2,  NP 
MA[  1,  ]  for  L  “  1  ,  NI 
T1 

(DX[L,  .1]  for  .1=1,  N)  for  L  ~  l ,  NI 
U)  ,  V  )  for  K  “  1 ,  NS 

II  for  11,  N 1 

Ul,  Ml,  (V  for  K  I  ,  Ml))  for  .1  1  ,  N2 

III,  V)  for  .11,  NP 
l  fur  I  "  1 ,  NR 


TABLE  2 


LIST  OF  VARIABLES  USED  IN  THE  INPUT  FUNCTIONS 


T,  corresponds  to  t 

X[I]  (1  =  1,  N),  corresponds  to  X^(t) 

FX[I,  Jj  (J  =  l,  Ml  for  each  1  =  1,  N),  corresponds  to 
f.  .(X) 

RT[I]  (J  =  1,  N),  corresponds  to  r^t) 


In  contrast  to  the  input  parameters,  which  are  entered  as 
numerical  values,  the  input  functions  f. . (X)  and  r^(t)  are 
entered  in  the  form  of  program  statements  in  the  procedure  FUNC. 
How  this  is  done  is  explained  in  Chap.  1. 


2 . 5  Limitations,  Timing 

Because  of  the  limited  capacity  of  the  computer  store,  the  upper 
hound  for  the  number  of  differential  equations  the  program  can 
solve  is  19.  However  the  practical  bound  is  lower  because  the 
running  time  for  19  equations  is  probably  excessive.  The  limitations 
resulting  from  the  output  format  ares 

Maximum  number  of  equations:  14 

Maximum  initial  value  of  X^s  9999  units 

Maximum  length  of  the  region  of  integration:  200  days. 

These  can  be  changed  if  necessary  by  altering  the  program. 

It  is  very  difficult  to  estimate  the  running  time  for  a  given  sot 
of  equations.  The  running  time  depends  on  the  number  of  equations, 
the  desired  bound  for  relative  accuracy,  the  form  of  the  equations, 
the  length  of  the  region  of  integration  and  the  output  form.  However, 
the  following  results  give  some  indication  of  the  running  time. 
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For  the  two  sample  problems  in  Chapter  3  the  number  of 
equations  was  2,  the  desired  bound  for  relative  accuracy  was  0.01, 
and  the  length  of  the  region  of  integration  was  100  days.  The  first 
has  one  contribution  to  each  attrition  function  f ^ (X) ,  and  the 
second  has  two  contributions.  The  running  times  were  5  and  11  minutes 
respectively . 

A  system  consisting  of  5  equations  and  a  maximum  of  two 
contributions  to  an  attrition  function  has  also  been  run.  The 
running-time  was  approximately  30  min. 
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3. 


SAMPLE  PROBLEMS 


Two  sample  problems  have  been  solved  to  illustrate  the  possibilities 
of  the  program.  The  first  one  shows  how  to  make  instantaneous 
introductions  of  units  and  how  to  change  the  analytical  form  of  the 
introduction  rates.  The  second  one  shows  how  to  manipulate  a 
system  with  more  than  one  contribution  to  the  attrition  of  some 
classes  of  units. 


3 . 1  Problem  1;  Lanchester  Linear  Law  with  Replacement- 

Consider  the  system  of  two  simultaneous  differential  equations 

Xj  «  -dX^X.  +  r, 

Xa  =  -0XjXa  +  r2 

with  initial  values  X*  (0)  =  1000,  Xa(0)=50j  and  attrition 

coefficients  a“5*l0“4,  0  =  2*lo"“^.  This  is  a  system  of 

Lanchester  equations  with 

fi  (x)  “  f u  (x)  =  ax^s 
f a  (X)  =  f'w(X)  =  0XxXs 

[only  one  contribution  f^.(X)  to  each  f \  (X)  J. 

An  additional  400  units  of  class  1  are  introduced  at  time  t  20. 
Introduction  rates: 


(t )  "  0  for  ail  t 


ra  ft) 


0  for  t  £  50 
0.1  for  t  >  50 


le 


The  following  part  of  the  optional  output  is  desired: 


a.  X2  Vs  Xj  (t  as  a  parameter) 


c.  ra  vs  t 


Input  parameters.  Let  the  region  of  integration  be  100  days. 
According  to  the  instantaneous  introduction  and  the  change  of 
introduction  rate,  the  region  of  integration  is  divided  into 
three  sub-regions:  0-20,  20-50,  and  50-100.  The  instantaneous 
introductions  must  be  specified  at  the  start  point  of  each  sub- 
region,  This  means  that  the  following  six  numbers  must  be  given 
as  input: 

Class  1  Class  2 

t  -  0  1000  .  50 

t  =--20  400  0 

t  "50  0  U 


The  time  interval  Dt  must  be  chosen  such  that  the  three  sub-regions 
become  an  integrsr  multiple  of  Dt .  Let  Dt  -  5.  This  implies  that, 
the  number  of  time  increments  Dt  in  each  sub-region  is  4*  b  and  10 
respectively.  Furthermore,  let  the  initial  value  for  the  step 
length  be  1  day,  the  lower  bound  for  the  step  length  be  0.01,  and 
the  desired  bound  for  the  relative  accuracy  be  0.01  (1$).  This 
gives  the  following  scheme  of  input  data  to  be  punched  on  paper 
tape  (for  definition  of  input  parameters  see  Table  L): 


DATA : 

2  5  1  0.01  0.01  1512100 

4  0  10 

0 

loon  50 
4(io  o 
o  i) 

2  1 

1  2 

O 


LND  OF  DATA; 


Input  functions.  The  input  functions  f^(X)  and  r^(t)  are 
easily  inserted  into  the  procedure  FUNC  by  copying  the  existing 
paper  tape  for  the  procedure  up  to  and  including  the  line 

''rnmniPnt  T n  qp  r't  thp  funrf  i  nnc  j  " 

and  punching  the  following  four  statements: 


FX[1, 1] :=0.0005*X[1]*X[2]  ; 

FX [  2 , 1  j  :  =0 . 0000  2*X[  1  ]*X[  2  ]  ; 

RT[1]:=0.0; 

if  L=1  or  L=2  then  RT[2]t=0.0  else  RT[2]j=0.1; 


(The  three  sib-regions  are  identified  by  the  values  Ls  1,  L-2 
and  L-  3  respectively.)  Then  copy  the  rest  of  the  tape.  The 
procedure  FUNC  then  takes  the  following  form: 


£Numerical ' solution  of  a  system  of  Lanchester  differential 
equations;?  ; 

begin  procedure  FUNC ( N , MI , T, X, F : FXS, FX, RT , L) ; 
value  N. MI; real  T {integer  N,MI,L; 
array  X,F,FXS,FX,RT; 
begin  integer  J,K; 

comment  Insert  the  input  functions; 


FX 

FX 

RT 


1,1 

2,1 


1] s  =0 .0  ; 


s=0.0005*X[l]*X[2] ; 
:=0.0000  2-”-X[l]-”-X[2j; 


if  L=1  or  L=2  then  RT[2]:=0 


0  el se  RT[2]:=0.1; 


for  J:=l  step  1  until  N  do 
begin  FXSf J] : =0 . 0 ; 

for  K:=l  step  1  until  MI  do 
FXSf J ] sFXSTjj+FXr J,Kl ; 
F[J]:=-FXS[JJ+RT[J]; 
end ; 

end  FUNC; 


Output .  The  numerical  values  of  all  the  output  variables  are  given 
in  printed  tables  obtained  from  the  output  paper  tape.  The  plotted 
output  is  presented  in  Figs.  2  to  5 • 
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FIG.  2  NUMBER  OF  UNITS  v*  TIME 


FIG.  3  NUMBERS  OF  UNIT  1  v*  NUMBERS  OF  UNIT  2 


Figure  2  shows  the  number  of  units  in  class  1  and  2  as  a  function 
of  time.  The  effect  of  the  instantaneous  introduction  of  units  in 
class  1  at  t—  20,  and  the  change  to  a  non-zero  introduction  rate 
in  class  2  at  t = 50  is  obvious.  However,  the  result  of  the  combat, 

if  continued,  does  not  he  romp  rloar  Pnom  fKj  »  c  i  - -  fwn  ■»  .  .  . 

-  -  *  me  pj.octJ.ng 

of  Xa  Vs  X*  with  time  as  a  parameter  in  Fig.  3  is  very  useful  in 
examining  this.  The  instantaneous  introduction  changes  an  apparent 
victory  for  class  2  to  an  apparent  victory  for  class  1.  However, 
after  changing  the  introduction  rate  of  class  2  from  zero  to  a 
constant  non-zero  value,  it  seems  as  if  class  2  will  win  the  combat. 

Figure  4  shows  the  number  of  units  lost  per  day  as  a  function  of 
time  for  class  I.  and  2.  The  two  curves  are  identical  in  this  case 
because  f;  is  proportional  to  f8  and  the  plotted  values  are 
given  as  a  percentage  of  the  initial  values. 

Figure  5  shows  the  introduction  rate  of  class  2  as  a  function 
of  time. 

3-2  Problem  2:  Lanchester  Mixed  Law 

Consider  the  system  of  two  simultaneous  Lanchester  differential 
equations 

*1  =  -(ttj  X:  X,  +as  xa)  +  r> 
xa  ~  -  ( Pj  Xx  Xa  +  3a  Xj, )  +  r2 

with  initial  values  X* (0)  -  1000,  Xa(0)  -  50,  and  attrition 
coefficients  **-4.10-4,  *,=0.1,  fj-l.S.lO"5,  |8  -0.002. 

This  system  has  two  contributions,  ^.(Xj.to  each  i\  (X) : 

‘VX)  -  fu(X)  +  f,?(X) 

fafx)  -=  r81ot)  +  ‘'BS(X) 
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where 


fjj  (x)  “  a,  x,  x, ,  fa(x)-«axt, 
fa(x)  =  6,  Xs  x8,  fM(x)  =  |8  xa. 

Let  the  introduction  rate  r} ft)  be  zero  for  all  t. 

The  following  part  of  the  optional  output  is  desired: 

a.  X2  vs  Xt  (t  as  a  parameter) 

b.  f .  vs  t,  i  =  1,  2. 

C*  ^  i  j  /  f  i  vst,  j  =  1,  2  for  each  i  =  1,  2. 

d.  fj/f,  vs  t. 

Input  parameters.  Let  the  length  of  the  region  of  integration,  the 
interval  Dt ,  the  initial  value  for  the  step  length,  the  lower  bound 
for  the  stop  length,  and  the  desired  bound  for  the  relative  accuracy 
be  the  same  as  in  Problem  1.  Since  neither  instantaneous 
introductions  nor  analytical  changes  of  the  functions  i\.(X)  or 
r^(t)  occur,  there  will  be  no  division  of  the  region  of  integration. 
This  leads  to  the  following  scheme  of  input  data  to  be  punched 
on  paper  tape  (for  definition  of  input  parameters  see  Table  1): 

DATA: 

2  5  l  0.01  0.01  2112021 

20 
0 

1000  50 
2  1 
1  2 

12  12 
2  2  12 
l  2 

END  OF  DATA; 
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Input  functions.  In  the  same  ua;'  a.?  in  Problem  1,  the  existing  paper 
tape  tor-  the  procedure  FI  NC  is  copied  and  the  following  statements 
inserted : 

F  a  ~  j .  j.  ! :  =-0.0004 :  xL  1  !  X[  -1 ;  FX[1,  2  i :  =0 ,1-xT  2  1 ; 

FXf  2,  l]  5=0.000018  Xf  L 1  * X l“  2]  j  FX[  2,  2] ;  =0 .002  >X[  2  1 ; 

RTF  1 1 r “RTf  21 ; -0 ,0 ; 


Output .  The  numerical  values  of  all  the  output  variables  are  given 
in  printed  tables  obtained  from  the  output  paper  tape.  The  plotted 
output  is  presented  in  Figs.  0  to  11. 

Figure  b  shows  the  number  of  units  in  class  1  and  2  as  a  function 
of  time. 

Figure  7  contains  the  plotting  of  Xz  Vs  Xx  with  time  as  a 
parameter.  This  figure  clearly  shows  the  result  of  the  combat, 
if  continued. 

Figure  8  shows  the  number  of  units  lost  per  day  as  a  function  of 
time  for  class  1  and  2. 

Figures  9  and  10  are  very  convenient  in  examining  the  relative 
importance  of  the  different  contributions  to  an  attrition 
function  1 (X) .  Figure  9  shows  that  fn  (  X)  constitutes  the 
greatest  part,  of  f  (Xi)  in  the  beginning  of  tin*  combat.  However', 
as  the  t,  i  me  increases,  1 1;.,  (X)  increases  while  1  U  (X)  decreases; 
consequent Ly  fu (X)  will  dominate  after  some  time.  The  same 
relation  appears  in  Fig.  lo  for  f2i  ( XI  and  • 

Figure  II  indicate^  that  number  of  class  1  units  lost-  per  class  2 
unit  lost,  increases  with  time,  emphasizing  the  fact  that  class  2 
will  win  the  combat  . 
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FIG.  10  PERCENTAGES  OF  LOSS  RATE  FOR  UNIT  2 
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FIG.  11  EXCHANGE  RATE 


General  Remarks  About  the  Input  Function 


Sample  problem  1  shows  how  to  change  the  analytical  form  of  one 

function  iv  ( t )  for  some  t.  The  same  is  of  course  possible  for 

anv  i'  l  t )  and  anv  f.  .IX).  If  some  change,  exist  s  for  a  func  tion, 
i  •  ij  — 

say  the  analytical  form  of  this  function  must  be  specif ied 

fur  every  sub-region  ( L  -  1,  N 1 )  in  a  statement  of  the  form 
ii  dil  uted  in  Problem  1. 


Since  MJ  is  the  maximum  number  of  contributions  1\^(X)  to 
f,(X)  ( i  *  1  n) ,  some  of  the  funct  ions  F  X  [1,  J  i  may  of  course  be 
equal  to  zero.  However,  all  of  them  have  to  be  inserted  into 
fUNCIJ}'  1.  MI  for  eacii  1  -1,  N).  Likewise,  all  the  functions 
Rf[  L  ]  ( l  -  1 ,  N)  must  be  inserted. 

If  it  happens  that,  two  or  more  of  the  contributions  f.  .(X.) 
tor  the  same  index  i  are  not  desired  as  output,  these  may  be 
put  together  in  one  function  FX[I,  J ] ,  This  will  reduce  the  runnin 
t  i me  on  the  computer. 


i 

i 
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APPENDIX  A 


C'OMPl'TLR  PRO(;k.AM  LISTINGS 


CNu— rleal  solution  of  a  system  of  Lanchester  dlffarantlal  equations?; 

begin  procedure  FUWC?N ,MI ,T  ,X ,P ,PXS ,FX  ,RT  ,LJ ; 

co— ant  This  procedure  contains  the  input  functions. 

Para— ters^ 

N,  number  of  aquations. 

MI,  maximum  number  of  contributions  to  ona  attrition  function, 

T,  a  variable  for  the  ti— , 

X,  array  containing  the  number  of  units  in  each  class  at  ti—  T. 

P,  array  containing  the  values  of  the  derivatives  at  ti—  T. 

PXS, array  containing  the  values  of  the  attrition  functions  at  ti—  T. 

PX ,  array  containing  the  values  of  the  different  contributions  to 

the  attrition  functions  at  ti—  T. 

RT,  array  containing  the  values  of  the  introduction  rate  functions 
at  ti—  T. 

L,  variable  which  has  the  value  1  when  T  in  first  sub-region,  2 
when  T  in  second  sub-region  ....... ; 

value  N  ,MI ireal  T; Integer  N  ,MI ,L; 
array  X ,F ,FXS  ,FX ,RT; 
begin  integer  J,K; 

com— nt  Insert  the  input  functions; 


for  Jr=1  step  1  until  N  do 
begin  PXSf jT:=0.0; 

for  K:rl  step  *  until  MI  do 
PXS[ Jl :=FXSrJl+PXr J , K ] ; 
FrJl:=-FXS[J]+RT[Jl; 
end; 

end  FUNC; 


procedure  MERS£»(N,MI  ,T,Z;H  HJKIN,E,X  ,F  ^S.FX.RT.LL^  ; 

coawnt  Tnia  procedure  perfora*  the  integration  fro*  T  to  2. 
Parameters  * 

Z,  end  point  ol  integration. 

H»  starting  value  for  the  atop  length. 

HMIN, lower  bound  for  tha  a tap  length. 

Ev  daalrad  bound  for  relative  accuracy,; 
value  W.M1 , HMIN ,E; 
real  T.Z.H.HMIV  .E- 
Integer  N.MI.LL; 
array  X ,? ,?XS .EX ,RT; 
comment  Variables: 

Q,  1.  fne  raat  of  tha  interval  of  integration. 

2.  ft  variable  used  in  tha  computation  of  the  fiva  aub-atapa. 

3.  5*?Truncatlon  error). 

'J.  1.  A  variable  uaed  In  tha  coaputatlona  of  the  fiva  aub-atapa, 

2.  iJiabaCXtl))  . 

H3=H/3 

riS ,  atcre#  H  for  a  poaalble  restoring, 

TS  ,  store#  T  for  a  poaalble  repeated  integration, 

ES-5*E 

I.J.aW.indax  variable*. 

G,  array  storing  X  for  a  poaalble  repeated  Integration. 
P„L,urraya  for  temporary  storing  In  the  coaputatlona  of  the  five 
aub-atepa. 

Boolean  variable#" 

BCvtrue  if  FO 

BF-falae  if  the  condition#  for  doubling  the  atep  length  are  not 
aatiaf led. 

BH=fal#e  if  the  atep  length  la  halved. 

Hfcrfalae  If  H  it  equal  to  the  reat  of  the  interval. 

BX-fala«  if  H=HM1W; 
begln  real  Q.U ,H3 ,HS ,TS ,E5; 

Intager 

array  G  ,P,Lh  :Nj; 

Boolean  BC , BE , BHBR ,BX ; 
switch  S SS : =SAVE , BACK ; 

comment  Check  eoae  input  parameters  and  Initialize; 

JUT  HMTN<0  then  HMIN:=O.OT*abs(H)  ; 

BH:^BR"=BX:=true; 

BC • -E<1  ; 

E:=absfE>  ; 

E5:=5*E; 

Hi=*bsOt) ; 

SAVE:  If  BC  then 

begir,  TS  *  -T ; 

for  J:=1  atep  1  until  N  do 

ofitl:=Xf  jf; 


BAcX:  H6:=H; 

Q:=TfH-Z; 

BE : =tru* : 

co— nt  Test  snd  of  integration  rang*; 

I-*  QX)  than 
begin  h:=Z-t; 

BH:=false: 

and 

Next  integrate  one  step; 

H3:=H/3; 

for  SW:=t  step  1  until  9  do 

begin  switch  SS< =SW1 ,SWEX ,SW3 ,SW4 ,SWAX ; 

FtJNC(N,W  ,T,X ,P ,FXS  ,FX  ,RT,1>L)  ; 
for  Is=1  step  1  until  N  do 
begin  switch  S :=ST1 ,ST2 ,ST3 ,ST4 ,ST5; 
switch  SSSS:=NEXT; 

Q:=H3»FCl1 rgoto  S[SW]; 

ST1:  L[I]:=U:=qjioto  NEXT; 

ST2:  U:=0.9*C<HLm>’;52i2  NBtTJ 
ST3 :  P[I]:=U:=3*Q; 

0:=0.37S*CU+Lm>  jgoto  NEXT; 

ST4:  L[I J :=U:=Lf I]+4*Q; 

U:=1 .3*(U-Pri]> ;goto  NEXT; 

STS:  U:=O.S*CQfLm>  ; 

Q:=abs(2*U-1 .5*(QfP[I]>) ; 

NEXT:  X[I] :=C[ I]+U; 
if  SW=S  then 
begin  if  BC  then 

begin  U:=abs(X[I]) ; 

U:=lf  U<»-3  then  ES  else  E5*U; 
cowant  Test  adjustnent  of  the  step; 
if  q>U  and  EX  then 
begin  BR:=true; 

BH:=false; 

H:=0.8*H; 
if  IKHMIN  then 
begin  H:=HMIN; 

BX:=falae; 

end  The  step  was  halved,  restore  X  and  T,  and  go  back 
for  repeated  integration  with  this  new  step; 
for  J:=l  step  1  until  N  do 

xfj]:=orjTV 

T:=TS jgoto  BACK; 
end  if  Q; 

if  q>0,03125*U  then  BE:=false; 
end  If  BC; 
if  not  BX  then 

begin  if  1=1  then  print  punch(l)  ,EEle10??; 

print  punch(l) .saaeline ,C*7 ,scaled(4) ,E*Q/U; 
end; 

end  if  SW; 
end  for  I; 
goto  SSfSirl; 
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SW1 :  T:=T+H3;goto  SWSX; 

SV3 :  1i*;-T»O.Sl*H3:|oto  S¥ZX; 

SW4 :  T:=7*0.3*B; 

SWAX:  H:=M;  comment  Dummy  stateaant; 

SID :  and  for  SW; 

4  #  »k.. 

begin  co—  nt  Test  a  possible  doubling  of  tha  atop; 
if  Bl  and  M  and  BR  than 
bagln  H:=2*Hi 
BX:=truo; 
and: 

BH:=trua; 
and  If  BC; 

if  BR  than  goto  SAVE: 

H:=HS; 
and  HER SON; 


procadura  AXPLTTi TM ,XM ,R1 ,R2,I1  ,12) : 

co—ant  Thia  procadura  plots  axas  and  axa  valuaa  for  tha  aaxlanm 
axa  valuaa  givan. 

Parameters : 

TV ,  maximum  absolsaa  value. 

XV,  maximum  ordinate  valua. 

Rl  ,R2, scaling  factors. 

II, 12, number  of  divltlona  raquirad  on  tha  tiro  axaa; 
real  TM.XM.RI ,R2; 

Integer  II  ,12; 
begin  intagar  J? 

II  :=rnif4.0)/10.0; 

R1 :=780. 0/<10.0*I1) ; 

I2s=rXMf4.0)/10.0; 

R2:=«00. 0/C10. 0*12) ; 

aatoriginC 200 ,R1  ,R2,1) ; 

axeaC10,10,I1 ,0,12,0) ; 

plottarCIO  ,1)  i 

for  J:=0  step  1  until  II  do 

bagln  movepenCJ*10„0-30.0/R1  , -20.0/R2)  ; 

print  dlgitsC3) ,J*10; 
and; 

for  J:=0  step  1  until  12  do 

bagln  aovapanC -90,0/RI , J*10.0-6.0/R2) ; 

print  diglts(3) ,J*10; 
end; 

and  AXPLOT; 


3  2 


Cl 


ommm nt  The  Min  prograa  contain*  *  loop  calling  MERSCN  npMttdlj  and 
storing  tha  values  of  tha  variables  at  the  and  points  of  ovary 
interval  DT.  Tha  rest  of  tha  Min  prograa  provides  for  output. 

Variables : 

DT,  ths  interval  length, 

TM,  end  point  of  tha  region  of  integration. 

XM,  YM.Mxiaua  axe  values  for  plotting. 
t  j  *  l  index  variables. 

M,  index  variable  for  the  ties  output  points. 

Ml ,N2,U,V, integer  variables  with  different  applications. 

Nl,  nuaber  of  sub-regions  of  integration. 

NS,N1 ,NR ,N2 .NP.inputparaietere  which  provide  for  the  optional  output, 

MA,  array  containing  the  nuaber  of  DT  in  each  sub-region. 

MK,  array  containing  the  value  of  M  at  the  start  point  and  end 
point  of  every  sub-region. 

NN,N3 ,N4,NM,N5,N6 .arrays  storing  the  different  Indies*  used  in  the 
optional  output, 

TPL,XPL,FXSPL,FXPL,RTPL, arrays  storing  the  values  of  T,X,FXS,FX  and 
RT  respectively  in  the  output  points. 

DX,  array  containing  the  instantaneous  introductions  of  units  at  the 
start  point  of  each  sub-region; 
begin  real  DT.H.HMIN.E.T.TH.Z.XM.TM.RI  ,R2; 

integer  I , J,K,L,M,N, II ,12,111 ,M1 .M2 ,NI ,NS ,H1 ,NR  ,N2  ,NP,U  ,V; 
read  N  ,DT  ,H  ,HMIN  ,E  ,MI  ,)f I  ,KS  ,N1  ,NR  N2,NP; 

J^f  N>  NF  then  I:=N  else  I:=NP; 

begin  Integer  array  MA[1 :NI],lTltl :2*NI] ,NNfC:Nl 1 ,N3 ,N4fO :H2] ,NMfO :N2  ,1 :MI1 , 
N5,N6[0:NPl; 

MisO; 

for  L:  =1  step  1  until  NI  do  begin  read  MA[L];M:=MtMA[L]  end; 

Ms=MfNI; 
punch(l)  ; 

print  ££1*1?NUMBER  OF  EQUATIONS ,  N  =? ,saaellne ,digit*C3)  ,N , 

£fil*1?A  GUESS  FCR  THE  STEP  SIZB.H  =? .saseline ,soaled(4) ,H , 

££lsl ? LONER  BOUND  FCR  THE  STEP  SIZE, RUIN  s? , sane line ,RMIH , 

££1*1 ?  DESIRED  BOUND  FOR  RELATIVE  ACCURACT,  E  =?  .saaeline  ,B; 
begin  array  X.F.FXS.RTCl :N] ,TPL[OjM]; 

ooassent  CBS  array  XPL[1 tM.l ill .FXSPL.RTPLfl :M,1 ;Nl .FXTl »M,1 :MI1 , 
FXPLn  !M,1  :N,1  JMI],0X[1  :NI,1  sN]; 
switch  S : =L1  ,L2,L3,L4,L5,L8,L7,L8,L9; 
eoswnt  Read  and  print  initial  values; 
read  T;TPL[0]:=T; 
for  L : -1  step  l  until  NI  do 
Tor  J:=1  step  1  until  N  do  read  DXTL.Jl; 

print  ££13s1? INSTANTANEOUS  INTRODUCTION  OF  UNITS:? ,££12s9?T?; 
for  Ital  step  1  until  N  do  print  saMllne  ,££s6?DX?  ,diglts(2)  ,  I ; 
TM:=T; 

for  L:=1  step  1  until  NI  do 

begin  print  EE1?? , see* line ,scaled(4) ,TM; 

for  J: =1  step  1  until  N  do  print  saaelln* .scaledO)  ,DXfL,J]; 

TM :  rTRkluftl  *DT ; 
end; 


print  C213e9?T?; 

for  I:=1  atop  1  until  K  do 

print  aaneline ,KBeT?X? .digitaf 2) , I ; 

for  J:=1  atop  1  until  N  do  begin  XPLCl ,J]:=SK[1 , J] ;X[ J] :=0.0  ond; 
M2:=Ot 

oo— nt  Start  integration  of  each  aub-raglon  with  now  initial 
valuaa  X[J); 

tor  L ; =1  otop  i  until  H 1  do 

begin  Ml  ;  =1(2+1  ;M2:=M1+MA[Ll ; MfT 2*1^— 1  ]:=M1  ;1M[2*L]  :=M2; 

TPLCM1 1 : =TPLf Ml  -1 ] ;print  CC17?  .aamellne ,aoaled(4) ,T; 
for  J : =1  atop  1  until  N  do 

bagln  Xf JT:=Xf Jl+DXfL , J] ; print  aaaoline ,acaled(5) ,XfJ]; 

.if  L>1  than  XPL[M1  ,  J]  :=100.0*XfJl/XPL[1  ,J]; 
end: 

oownant  Call  FCNC  for  computation  of  the  output  variablaa 
at  the  atart  point  of  aub-region  no.  L; 
rWCCM.Ml  .T.X.F.rXS  FX.RT.L); 
for  J : =1  atop  1  until  N  do 
bagln  FXSPLtm , jl :=PXS[j77RTPLrM1 ,J]:=RT[J]j 
for  K:=1  atop  1  until  MI  do 
FXPL[M1 , J,K]: =rX[j  ,K] ; 
and: 

ooaaaant  TM-the  ralua  of  T  at  the  and  point  of  aub-raglon  no.  L; 
TM:=T+MA[Ll*im0.1  ;M:=M1 ; 

coament  Call  MBtSON  and  PUMC  repeatedly  for  computation  of  the 
output  variablaa  in  each  output  point  in  aub*reglon  no.  L; 
for  Z : -T+DT  atop  DT  until  TM  do 
begin  MXRSOMCN ,MI  ,T ,Z  ,H,HMIN,E,X,r,rXS, FX.RT.L)  ; 

“tCrcc  r  ,mi  .  t  ,x  ,r ,  ns  ,rx ,  rt  ,l>  ; 

M:=M+1 ;TPL[Ml:=T; 

comment  Print  T  and  X  in  output  point  no.  M.  Store  the  name 
value*  for  l*tar  plotting; 

Print  efil??.aaneline,acal«d<4> ,T; 

for  J : =1  atop  1  until  R  do 

bagln  print  aamellne  .acaledO) ,X[J]; 

XPL[M,J] :=100.0*X[ Jl/XPL[1 ,Jl; 
PXSPL[M.Jl:=FXS[Jl;RTPL[M,J]:=RTrJl5 
for  K : =1  atop  1  until  MI  do 
FXPUM,  J,K]:=FXrJ,Kl; 

•"d: 

end: 

end  for  L; 

co— nt  Plot  X[  1 1 , 1=1  ,K ,  va  tine; 

XM:  =100.0; 

for  J:=1  atep  1  until  N  do 
for  L:=2  atep  1  until  2* NI  do 

Tggln  Ml  :'35rL];if  XPLfMI  ,jT>XM  than  XM:=XPI.[M1  .Jl;ond; 
AXPLOTfTM.XM.R1 ,R2,I1 ,12) ; 
movepen(390.0/R1 .-40.0/R2); 
print  plotterflO.1) .CTiaefdaye)?; 
aervepenf  -50 .0/R1  .100.0/R2) ; 

print  plotterOO  ,3)  .BHuaber  of  unite  in  par  oent  of  lnitlalvalue? ; 
move pen ( -1 90.0/R1  .900.0/R2)  ; 

print  plotterdO  ,1)  ,einitial?;aovepen(-180.0/R1  .480.0/R2); 
print  plotterOC ,1)  .Cvaluee:?; 
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for  J:=1  «t»p  1  until  N  do 

bagin  aovapanf -180.0/M ,(480.0-20.0*  J)/R2)  ; 

print  plottardO  ,1)  ,EX?  ,digits(  ?)  . J ,£=? ,sl ignad(4 ,0)  ,XPUl  ,J]  ; 
and; 

for  J:=1  stop  1  until  N  do 
bagin  ■ovapancri'Li  1  J , 1 UO .  o) ; 
panlowar; 

for  I: =2  itip  1  until  M  do 
dmwlin#(TPLTl]  ,XPLn  ,Jl)  ; 
print  plottarCIO  ,1)  ,EX?,digits(2)  ,J; 

•nd  of  plot  X[I] ; 

■ovapan(0 , 20* 12) ; 

co— nt  Plot  X[I]  vs  xm  for  spaclflad  lndax  pairs  (I  ,J)  ; 

Tf  NSaO  than  goto  LI j 
for  K:=1  atap  1  until  NS  do 
TSagin  road  U,V; 
for  Ml : =V  do 
bagin  for  M2:=U  do 

bagin  XM:=TM:=100.0; 

for  L:=2  it op  1  until  2*NI  do 

bagin  if  XPLL MMf L  ]  .  M  ]  >XM  than  XM:=XPL[l«irL)  ,M1  ] ; 

If  mrnUfLl  , M2 ]>YM  £han''TB:sXPLtMM[L]  ,1(21; 

•nd; 

AXPL0T(XM ,  YM ,  R1  ,R2,I1  ,12) j 
aovapan(100.0/R1  .-40.0/R2)  ; 
print  plottardO  ,1)  , 

ENuabar  of  units  in  p«r  cant  of  initialvalua ,  X? ,digits( 2)  ,M1 ; 


aovapan(-30.0/R1 .30.0/R2) ; 
print  plottardO. 3)  , 

ENuabar  of  units  in  par  c#nt  of  initialvalu# ,  X?,dlglts(2) ,M2; 
aovapan(-1 80.0/R1 .900.0/R2) ; 

print  plottardO  ,1) .Elnltisl? ;aovapan(-l 80.0/R1 .480.0/R2) ; 
print  plottardO, 1) ,£valuas:?; 
for  J : =1  ,2  do 

bsgln  aovapan(-1 80.0/R1 ,(480, 0-20, 0*J)/R2) ; 

print  plottar(10,1)  ,£X? ,digits(2)  ,(i_f  J=1  than  Ml  also  M2),£=?, 
alignad(4,0)  ,XPLtl  ,(if  J=1  than  Ml  also  M2)]; 

and; 

■ova pan  (100 ,100)  ; 

print  plottardO, 1>  ,£  T=?  ,alignad(2,0)  .TPLH  1 ; 

aovapanClOO  ,100)  ; 

panlowar; 

for  I : =2  stap  1  until  M  do 
drawl lna(XPLf I ,M1 ] ,XPLf I  ,M2l) ; 
aovapan(XPLtM ,M1 1-60.0/M 1XPLrM,M2]> ; 
print  plottardO  ,1 )  ,£T=?  ,»1  lgnad(3 .0)  .TPLfM]  ; 
aovapan(0,'0*I2) ; 
and ; 


co—  nt  Print  FXSfl]  va  tiae  for  apoclfiod  lndlc**  I; 

punch(l)  ;J^f  NlrO  th*n  goto  L2; 

print  ££13a9?T? ; 

for  I:=1  *t*p  1  until  N1  do 

begin  road  Uj 

WWfH:=U;print  a— line ,££»7?F? ,digita(2) ,U; 

»na ; 

for  J:  =  1  *t*p  1  unt 1 1  M  do 

Bogin  print  ££122  .aoaolin* .«c»l*dC4> ,TPL[ J} ; 
for  I:  si  atop  1  until  N1  do 
print  aaaallne ,acolod(5) .PXSPLf J  ; 

ond  of  print  FXSril; 

co— ont  Print  FXf I , J]/7XS[ 1]  for  *p*cifi*d  lndico*  I,  ond  specified 
indlco*  J  for  ooch  I; 
if  K2=0  than  goto  L3; 
print  ££13*9?T? : 

Tor  J : =1  atog  1  until  N2  do 
Sogln  rood  U,M1 ;N3fJ] :=U;H4r J]:=M1 ; 
for  K : si  atop  1  until  Ml  do 
bogln  rood  V; 

Print  soaolin*  ,££a42F2 ,diglta(2) ,U ,V; 

NHrJ.K]:=V; 

for  I : =1  atop  1  until  M  do 

FXPL[I ,U.V] : =100.0 *FXPLf I ,V ,V]/TXSPh[ I  ,U] ; 

ond; 

print  ££1*1022; 
ond ; 

for  I ; =1  atop  i  until  M  do 

Bogin  print  ££122 .aaaolino ,acalod(4) ,TPL[IJ ; 
for  J : =1  atop  1  until  N2  do 
bogln  Ml 

for  K : si  atop  1  until  N4[J]  do 
bogln  M2  =KMU  , kTJ 

print  ao—  1  in*  .acolodfS)  ,FXPL[  I  ,M1  ,M2] ; 
ond; 

print  ££1*1022 ; 
ond; 

ond  of  print  FX[ I , J]/FXSr I] ; 

co— nt  Print  FXSni/FXSfJ]  for  apocifiod  index  polra  (I,J); 

if  HP=0  than  goto  L4, 

print  ££13a9?T?; 

for  J:=1  atop  1  until  VP  do 

bogln  road  U  ,V;N3[Jl : =U;N6lJl :=V; 

print  a— lino, ££*'’???  ,digita(?)  ,U,£/F?,V; 
for  I :  =1  atop  1  until  M  do 
XPU  I  ,  J]  :  sFXSPUT7uT/rXSPL[  I  ,V]  ; 
ond; 

for  1=1  atop  1  unti 1  M  do 

bogln  print  ££12? , aaaol ino ,acalod(4) ,TPLf I] ; 
for  J= I  atop  1  until  NP  do 
print  aajoel  ino  ,acolod(3)  ,XPL[  I  ,Jl ; 
ond  of  print  FXST I  l/PXSf  *11 ; 


cniraer.t  Plot  rxsm  vs  time; 

L4:  If  thi-i  got  r  1.5; 

XM :  r'00  : 

for  J 1  step  unt  i  1  N 1  do 
begin  IT!  :=KN[  J]  ; 

tor  I:-4  “ ct»p  j  until  m  do 

KXSPLr  i  .Ml  1 :  =100.0*FXSPlT7  ,M1  3/FXSPlT  1  .Ml-); 
for  I.:  =2  stop  1  until  2*NI  do 

O eg.m  m2  -MM[l.l;lf  FXSPL[ M2 , Ml ]>XM  then  XM : -FXSPLf  M2  ,M1 1  end; 
end; 

AXPLOT<TM,XM,RI  ,R2,IH  .12)  ; 
sovoper.(35C.O/Rl  .-4G.0/R2)  ; 
print  plctterOO  ,1)  ,£Ti«e<davs)  ; 

»ovspe:V-50.0/R1  ,'0.0/R2>  ; 

print  plot  ter  r  0 ,3)  .CN'umb .  of  units  lost  per  unit  time,  in  per  cent 
of  init.val.?; 

m&vepsn'C-j  80,0/Pi  ,500.0 /R2;  ; 

print  piotterOO  .')  ,EInitial?;s»ovepen(  -ISO-O/RI  .480.0/R2)  ; 
nrir.t  plot  ter  f  1  0  ,  •  '  ,£vftlues:?; 
for  J  :  =";  #  t «-p  1  until  N1  do 

begin  noveptin<  *1  80  ,0/Rl  ,  (480 .0-20 ,0*J) /R2)  ;M1  :=NN[J]  , 

print  plotter<l0,1>  ,EF?  , digits? ?>  Ml  .€=■?  ,aligned?2,r>  .FXSPLM  ,Mll 

s  nd ; 

Tor  J--1  *'_*.■•£  i  until  N1  do 

neg in  mov<?piv  i {  I'Pt.r  "*Y". '■  00 ,  oY  •  M1  ;  =NNf  Jl  ;  penlower ; 

for  ).  =2  stft£  untlj.  M  do  drawl  ino(TPLm  .FXSPLr  I  .Mil?  ; 
pTi r  t  pir/tfer; ,  <rT)  £F?, digits?  2)  ,M1  ; 
of  plot  rXafll; 
mOvopen-'O  .  20*  (2^  ; 
common t  Plc-t  FXl"  I ,  Jl/FXSf  I-] ; 

1.3:  if  K  2.0  then  goto  L6 , 

for  3  fop  1  until  N2  do 

begin  U:=N3f ul;M1 ;rN4f Jl;XM:=0.0; 
for  K:=1  stop  1  until  Ml  do 
begin  V=VMrJ,Kl; 

'or  !.:=1  step  1  unt  i  1  2*NI  do 
toegir  M2:=MM[l.l;H:=abs<FXPlJM21l)1Vl>  ; 

it  '-"'XM  then  XM  -H; 

_*nd  . 
end , 

AXf-ljOTT  1M  ,XM  ,R1  ,R2,T  ;  ,  T2>  ; 
m.ivopen :  J  V) .  0/R1  ,  -40 . 0  /ri2>  ; 
print  plotter?' 0  I'  .FTimeCdays) 
movepenf  -  50 .0  /R1  ,  '00.0/R2)  ; 

print  plotter? 10 ,3>  ,CF  I  J 'T  I  ,in  per  cent?; 

T_ or  K  :  -1  stop  1  unt  i  1  Ml  do 
begin  v'  -NMT.J  ,Kl’; 

JJ  Fxr :  r  ;  .11 .  /V0  then 

begin  ntovepeniTPI.r’ 1  , -FXPl.r  1  ,U  ,Vl)  ;penlower : 

for  I  . -2  s t op  1  until  M  do  drawline<TPL[  II  ,-FXPI.r  I  ,U  ,Vl)  ; 
piOt  terdO  ,i )  ,f  -F?7<1igit  *C'’3  ,U  ,V  ,C  -T?  ,1’i 

end  else 


bagln  aovap#n'{TPL[l ] ,FXPL[l  ,1! ,V]) ; panicraar ; 

for  I:  =2  atap  1  until  M  do  drawlina(TPL[ I ] ,FXPL[ I , U ,V]) ; 
print  plottarOO,’)  ,  E-t-F?  ,digit«(2) ,U ,V .E/F? ,U; 
and ; 
and ; 

_•  ^  h  nM»vn^ 

w  »  aprtfii  \  u  ,  v#  *«•/  ( 

•nd  of  plot  PXn,J]/nCS[l1i 
co<— nt  Plot  rxsnJ/rxstJ]; 
if  MP=0  than  goto  L7 ; V . =0 ;N * =HP ; 
for  J :  =1  atap  1  until  N  do 
bagln  XM : =0.0 ; 

for  L i =1  atap  1  until  2*NI  do 
bagln  Ml  :=MMTM 

it  XPLtMl  , J 1  >XM  than  XM:=XPL[M1  ,Jl; 

•nd; 

E  =1 .0; 

U  XM> 200,0  th»n 
bag in  K; =1 ■ 

for  K-rK*10  whtla  XM>200.0  do 

bagin  M2  :=K;XM:=XM/10.0  and ;E ; =1 .0/M2; 

•  nd; 

"if  XMC20.0  than 
bagln  K ■ = 1 ; 

for  K  i  =K*1 0  ah  1  la  X*r20.0  do 
bagln  (H?  iK:XM:=XM*10.0  and  E  =M2; 

•  nd ; 

AXPLOT^TM.XM.R .  R2 , II ,12); 
oovapanf  350  .0/R1  ,  -40 .0/R2)  , 
print  plotterOO  ,1)  ,€Ti**<d»y«>?: 
nova pan( -50.0/R1  ,  200.0/R2)  ; 

_l_f  V=0  than  print  plott«r<10 ,3)  ,  f  r*«pointf  5)  ,E,E*r  I  /  F  J? 

•  !••  print  plottarOO .  3) , fr*apoint( 5) ,E ,£*R?  .digital 2) , NST J 1 ; 
■ovap#n(TPLh ] ,E*XPLfl  ,  J]) ;p«nlow»r; 

for  1-2  atap  i  until  K  do  drawlinafTPLtll  ,E*XPLf  I  ,J])  ; 

JLf  V=0  than  print  plottarCIO  ,1  >  ,£F? ,digit«(?)  ,N3[J],£  /  F?,N6[J] 
alaa  print  plott*rv'10  ,l)  ,£R?  ,diglta<2)  ,N5[  J]  ; 

*ovap*nfO,'X)«I2)  ; 
and  of  plot  FXST H/FXSf ll ; 
if  V=1  than  goto  L8; 

coaraant  Plot  RtT  1 1  vi  tl»a  for  apacifiad  indicaa  I; 
if  KR=0  than  goto  1,8 ; V  ;  =1  ; N  ’  rNR ; 
for  J : =1  atap  1  unt i 1  NR  do 
bagln  raad  U;N3[ jl :=U; 

for  I ; =1  atap  1  until  M  do  XPLf I .3] : =RTPL[ I ,U] ; 
and;  i^o  to  1,9; 

coMBant  and  of  plot  RTf  1 1 ; 

N  =N;  cowant  IXu#*y, 


